Shirley Interpolation - Optimal basis sets for detailed Brillouin zone integrations 

revisited 



David PrendergaslQ 

The Molecular Foundry, Lawrence Berkeley National Laboratory, Berkeley, California 94720 



o 
o 

(N 

Oh' 
(D 
00 



C/3 



Steven G. Louie 

Department of Physics, University of California, Berkeley, California 94720 and 
The Molecular Foundry, Lawrence Berkeley National Laboratory, Berkeley, California 94720 

We present a new implementation of the k-space interpolation scheme for electronic structure pre- 
sented by E. L. Shirley, Phys. Rev. B 54, 16464 (1996). The method permits the construction of a 
compact k-dependent Hamiltonian using a numerically optimal basis derived from a coarse-grained 
set of effective single-particle electronic structure calculations (based on density functional theory 
in this work). We provide some generalizations of the initial approach which reduce the number 
of required initial electronic structure calculations, enabling accurate interpolation over the entire 
Brillouin zone based on calculations at the zone-center only for large systems. We also generalize the 
representation of non-local Hamiltonians, leading to a more efiicient implementation which permits 
the use of both norm-conserving and ultrasoft pseudopotentials in the input calculations. Numer- 
ically interpolated electronic eigenvalues with accuracy that is within 0.01 eV can be produced at 
very little computational cost. Furthermore, accurate eigenfunctions - expressed in the optimal basis 
- provide easy access to useful matrix elements for simulating spectroscopy and we provide details 
for computing optical transition amplitudes. The approach is also applicable to other theoretical 
frameworks such as the Dyson equation for quasiparticle excitations or the Bethe-Salpeter equation 
for optical responses. 
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I. INTRODUCTION 

Providing efficient access to accurate electronic struc- 
ture is vital to accelerating research in materials science 
and condensed matter physics. This can be achieved di- 
rectly by increasing the availablility of computational re- 
sources, by developing faster numerical methods, or by 
switching to more compact numerical representations. 
However, if one adheres to existing methods and repre- 
sentations, one might reasonably ask if more information 
can be extracted efficiently from such approaches. In this 
work, we outline an efficient approach to extracting de- 
tailed information on electronic structure for arbitrary 
electron wave vector k. This method applies not only 
to periodic systems - where k is well-defined - but also 
to models of aperiodic systems within the supercell ap- 
proach For example, periodic calculations are often 
used to simulate isolated molecules in large supercells and 
disordered condensed phases are commonly modeled us- 
ing a supercell of sufficient size to contain relevant struc- 
tural features. Providing efficient access to first princi- 
ples electronic band structure and matrix elements over 
the entire Brillouin zone (BZ) supports a wide range of 
research topics from Fermi surface exploration in super- 
conducting materials to detailed simulated spectroscopy 
of dispersive bands or high energy excitations. 

In 1996, Shirley [2] outlined an approach within effec- 
tive single-particle electronic structure for constructing 
an optimal basis which spans the BZ and can be used to 
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build a compact k-dependent Hamiltonian based on some 
coarse-grained reference calculations. He applied this 
approach in detailed explorations of the dispersion and 
spectroscopy of crystalline systems: siHcon, germanium, 
graphite, hexagonal boron nitride, lithium fiuoride, and 
calcium fiuoride. The efficacy of his approach was tested 
by examination of electron band structures, densities of 
states, dielectric properties, x-ray resonance fiuorescence 
and incoherent emission spectra, and photoelectron spec- 
troscopy, with details provided or referenced in [5]. He 
also used this basis in developing efficient approaches to 
Bethe-Salpeter calculations for valence-band [3| and core- 
level spectroscopy To our knowledge, this advanta- 
geous approach has not been widely appHed outside of 
Shirley's research. In this work, we outHne a generalized 
implementation of Shirley's interpolation scheme, which 
has been incorporated as a post-processing tool for use 
with the Quantum-ESPRESSO ^ open-source electronic 
structure package. We hope that the outline provided 
here indicates how easily Shirley interpolation might be 
implemented in other electronic structure codes. Fur- 
thermore, we illustrate that the advantages of such an 
approach are clear when applied to large supercell cal- 
culations, where zone-center electronic structure calcula- 
tions are sufficient to accurately reproduce the electronic 
structure throughout the BZ. This particular implemen- 
tation has already been used in simulating x-ray absorp- 
tion spectra of molecules in vacuo @, 0|and in solution 
i- 

Another commonly used interpolation scheme for elec- 
tronic structure exploits maximally localized Wannier 
functions This approach builds a set of Wannier 
functions to describe a band complex and minimizes 
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their spread in real-space. The resulting functions 
can be quite localized and enable the calculation of 
first-principles tight-binding parameters for use in cal- 
culations of Berry's phase polarizationll d lllll- elec- 
tron transport (T2|, electron-phonon coupling[l3l.ll4l|- and 
much more. The Shirley interpolation method does not 
produce spatially localized functions. However, for inter- 
polation purposes, we show that it is much more auto- 
matic to use and computationally less expensive in terms 
of required input calculations. Also, conduction band 
states may be generated as easily as valence band states. 
In particular, there are no intrinsic difficulties in treating 
metallic systems and no specific requirements for disen- 
tanglement of dispersive bands [l5|. 

This paper is organized as follows: In SectionlUlwe pro- 
vide a summary of the work outlined in Ref . In Section 
Illll and lrVl we focus on the advances in our particular im- 
plementation over the original work. Section |V] focuses 
on some applications which highlight the advantages of 
the approach. Section IVTl differentiates Shirley interpola- 
tion from Wannier interpolation. Finally, in Section IVIII 
we provide some potential applications of this approach 
and then summarize our conclusions in Section IVIIII 



II. BACKGROUND TO THE METHOD 

A brief summary of Shirley's approach is provided here 
to establish the context of our own work. The process of 
building a compact k-dependent Hamiltonian in the op- 
timal basis for Brillouin zone sampling is outlined in Fig- 
ure[TJ We assume an effective single-particle Hamiltonian 
(based on Kohn-Sham density functional theoryfll, 
in this work). A self-consistent charge density is gener- 
ated by sufficient sampling of the Brillouin zone. Then, 
if necessary, a set of states is calculated from this den- 
sity for a user-specified set of band indices and k-points. 
The periodic parts of these Bloch-states are extracted 
and used to construct an overlap matrix which is then 
diagonalized to isolate linear dependence in this basis of 
periodic functions. By ordering the overlap eigenvalues 
by decreasing magnitude, we may select the optimal ba- 
sis subject to a user-defined tolerance. In our approach 
we truncate the basis by specifying a tolerance e for the 
neglected fraction of the trace of the overlap. 

Since we have removed the plane-wave envelope func- 
tions from the Bloch-states, a k-dependent Hamiltonian 
is required, and we represent this Hamiltonian H(k) in 
the optimal basis. Each component of the Hamiltonian is 
expanded as a polynomial in k. The kinetic energy has an 
analytic quadratic form, as indicated in Fig. [TJ Specifi- 
cally, with access to the Fourier coefficients of the optimal 
basis functions Bi{G), if one expands the k-dependent 
kinetic energy operator, one obtains: 
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Figure 1: The steps involved in building the compact k- 
dependent Hamiltonian in the optimal basis, beginning (top) 
with a self-consistent charge density, from which states are 
generated for a range of bands and k-points. An overlap 
matrix is constructed from the periodic parts of these Bloch 
states and diagonalized. The optimal basis is chosen as those 
eigenvectors of the overlap matrix (ordered by eigenvalue 
magnitude) which span a user-defined fraction of the space 
defined by the input states. The k-dependent Hamiltonian is 
constructed in this basis from its various parts: kinetic energy, 
local potential, and non-local potential. 
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The local potential is constant with respect to k and its 
matrix elements are efficiently computed within a plane- 
wave code using Fast Fourier Transforms. The non-local 
potential requires fitting on a grid in k-space. In Shirley's 
implementation [3], k-dependent matrix elements of the 
non-local potential operator were evaluated at points on 
a uniform Cartesian grid which contained the entire first 
Brillouin zone and these were fitted to polynomials in k 
to enable interpolation between the points. Explicit ex- 
pressions for each term in the Hamiltonian were provided 
in the original work. 

The outcome of these steps is a set of coefficients which 
can be used to construct the matrix i?(k) for any k and 
then diagonalize it to produce the eigenvectors and eigen- 
values in good agreement with an equivalent solution to 
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the underlying Hamiltonian at the same k-point. The 
advantage of Shirley's approach is that one reduces the 
size of the problem to be solved (i.e., the dimension of H) 
such that detailed explorations of the eigenspectrum be- 
come tractable. For instance, one could expect to switch 
from thousands of plane-waves to perhaps tens of opti- 
mal basis functions, thereby reducing a tough iterative 
(most likely parallelized) diagonalization to a trivial di- 
rect diagonalization which can be solved efficiently on a 
single processor. Shirley's satisfaction with the efficacy 
of this approach is clear in his original paper and also 
from the large number of detailed spectroscopic simula- 
tions enabled by it. 



III. MORE EFFICIENT ROUTES TO THE 
OPTIMAL BASIS 



A. Choice of k-points for the optimal basis 



The original scheme outlined particular uniform Carte- 
sian grids in k-space at which eigensolutions were gener- 
ated using the DFT code of choice and provided as input 
to construct the optimal basis by diagonalization of their 
overlap matrix. Henceforth, we shall refer to these DFT 
eigensolutions as "input states" and the k-space grid on 
which they are calculated as the "input grid". In gen- 
eral, the "input grid" need not coincide with that used in 
the original self-consistent DFT calculation, which gen- 
erated the self-consistent charge density, and the range 
of band indices for a given application may also differ, 
particularly when exploring the unoccupied spectrum. 
Therefore, the input states are often generated directly 
as solutions to the Kohn-Sham equations for a fixed self- 
consistent field. In the original work, the input grids 
contained the entire first Brillouin zone, with the aim of 
reproducing the eigensolutions at all points in that zone. 
It was not clear from the original work how the final ac- 
curacy would be affected by particular choices of such 
coarse input grids. Furthermore, the choice of input grid 
varied with lattice symmetry due to the differences in 
the shape of the first Brillouin zone in Cartesian space. 
In this work, we instead present a more general and au- 
tomatic approach to sampling k-space based on uniform 
grids in reciprocal lattice space, spanning the unit cube 
[0,1]^. This means that for any lattice symmetry, the 
input grid of k-points may be characterized uniquely by 
three integers ni x n2 x 713, much like a typical electronic 
structure calculation. Furthermore, one can be sure that 
this grid spans the volume of the Brillouin zone. In the 
original scheme, the use of a Cartesian input grid leads 
to the inclusion of some k-points lying outside the zone 
boundary for non-orthorhombic cells. 



B. Building in periodicity with respect to k 

As stated by Shirley, the k-dependent Hamiltonian 
does not impose periodicity in k-space, and so, one must 
be careful regarding the k-points one passes to the Hamil- 
tonian for diagonalization. We illustrate this point in Fig- 
ure [2] by examining the band structure of bcc Na along 
the A fine connecting the F and H points, with exten- 
sion to the 2H point, which is equivalent to F. We chose 
Na since its electronic structure is well described at the 
DFT level using a local pseudopotential, thereby remov- 
ing any complications associated with interpolating the 
non-local potential. In this calculation, we employed a 
norm-conserving pseudopotential and a 30 Ry kinetic en- 
ergy cut-off. The interpolated band structure is clearly 
dependent on the choice of input k-points used to gener- 
ate the basis. Using the F-point alone [Fig. Elja)] results 
in accurate electronic bands at that point, but large er- 
rors as one follows the A fine. Most notably, the periodic 
image of the zone-center, 2H, is completely wrong, which 
emphasizes that we have no explicit periodic boundary 
condition in our k-dependent Hamiltonian. Inclusion of 
the zone-boundary H point [Fig. EJb)] leads to marked 
improvement in accuracy along A but leads to some in- 
accuracy once we leave the zone-centered first Brillouin 
zone - again the 2H point is not reproduced. However, 
the ability to obtain excellent agreement in band struc- 
ture between the input grid points naturally prompts one 
to continue adding points to enable reproducibility over a 
larger region of k-space. Explicitly adding the 2H point 
[Fig. Etc)] does indeed almost restore the correct sym- 
metry of the band structure, albeit only in the neigh- 
borhood of this [F, 2H] interval - we should expect no 
reproducibility outside this interval. At this point, one 
can appreciate Shirley's original Brillouin zone-spanning 
choice of input grids, as they guarantee accurate repro- 
ducibility of band structure throughout the zone. How- 
ever, for cases where certain high-symmetry points are 
not included in the input grid, inaccuracies may appear. 

We have observed that the input grid does not define a 
"fit" in the usual sense of interpolation, with reproducibil- 
ity decreasing in accuracy as one explores points farther 
away (in the Cartesian sense) from the grid points. In 
fact, when we include the H point, we notice that the 
entire A line is accurately reproduced. Furthermore, if 
we were to include F and its periodic image 2H alone, 
we would see reasonable reproducibility across the en- 
tire A fine. Clearly, it seems that there is sufficient k- 
dependence built into H(k) to accurately describe the 
full-zone, once we provide additional constraints on the 
symmetry via the input states. This is equivalent to a 
phase-constraint for the optimal basis, as we shall see 
next section. And so, in our implementation we pro- 
vide input grids chosen uniformly from the unit cube, in- 
cluding all corners of the cube. We furthermore impose 
periodicity on k-points requested for diagonalization by 
mapping them first to the unit cube (in reciprocal lat- 
tice coordinates) , since we have no guarantee of accuracy 
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Figure 2: Accuracy in reproducing the band structure of bcc Na with respect to choice of input k-point grid. Band energies 
are reported in eV with respect to the Fermi level and k-points follow the A line from the zone center (F) to one of the zone 
boundaries {H) and beyond to a periodic image of the zone-center (277). DFT (Shirley interpolated) band structure is indicated 
by black lines (red dots). The size of the input k-point grid is varied: (a) F only; (b) F and H; (c) F, H and 2H; (d) F together 
with its seven periodic images from the corners of the unit cube [0, 1]^. 



outside the three-dimensional interval [0, 1]'^. 



C. Accurate interpolation using just the F-point 

If we choose our input grid from the unit cube, it seems 
clear now that we should include all corners of the cube 
- r and its seven periodic images in three dimensions. 
Always wanting to reduce computational eflFort, we im- 
mediately see that the input states originating from these 
periodic images differ only in phase from those at P. Fur- 
thermore, we will see that, for large supercells (small 
Brillouin zones), it is sufficient, using a DFT calculation, 
to generate input states at F only, generating periodic 
images of these states at the corners of the unit cube 
in a small number of steps. Figure [2^d) illustrates how 
well this approach works for bcc Na. The resulting band 
structure is indistinguishable (by eye) from the original. 
The root-mean-square error along the indicated path is 
5.5 meV. We may reduce this error by including more 
input k-points. However, we note that this error will not 
tend to zero, since it is comparable to the error associated 
with changing the kinetic energy cut-off and is related to 
a slight inconsistency in the number of plane-waves be- 
tween periodic functions from different k-points (as noted 
by Shirley). In our implementation, we retain all wave 
vectors G, such that i|k -|- G|^ < Ecut, padding with 
zeros the coefficients of those functions with missing G. 
Note that numerical differences related to this effect re- 
duce in magnitude for larger cut-offs, but in essence they 
are inconsequential, given that the original DFT calcu- 
lated eigenvalues would change as much upon varying 
Ecut- Generally, for large supercells we find that F-point 
sampling is adequate to reproduce all of the band struc- 
ture with an accuracy that is within 10 meV. 



To reduce the cost of the input DFT calculation, 
we construct the periodic images of the input states. 
The transformation of a given periodic function w„k to 
Unk+Go is easy to obtain in plane-wave representations. 
Since periodicity implies that 

e-**=^°'"u„k(r) = u„k+Go(i-) 
then, expanding in Fourier coefficients, we have 

^c„k(G)e'(«-°'')'- = ^c„k+Go(G)e^«-, 

G G 

which implies that the Fourier coefficients are 
ultimately reordered according to c„k(G + Gq) — 

Cnk+Go(G!-). 

For applications to large systems, where wave functions 
are necessarily distributed over many processors, such 
reordering may be complicated to implement. In this 
case, a simpler, albeit less efficient approach, is to exploit 
the native implementation of the Fourier transform to 
make such a transformation. Suppose that we start with 
the wave functions in Fourier space. Then we follow this 
map 

c„k(G) Unk{r) 

I 

c„k+Go(G) ^ e-*G°-'-u„k(r) 

where we first back-transform to real-space, then mul- 
tiply by the function e~''^° '' for each r, and then Fourier 
transform again to reciprocal space. 

Note that for cases where the F-point alone is insuf- 
ficient, this scheme could also be generalized to expand 
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input states to the star of a given input k-point by em- 
ploying the Httle-group of that k-point as determined by 
the lattice and atomic basis symmetry. 

IV. GENERALIZATION OF THE NON-LOCAL 
POTENTIAL 

A. Generalized Kleinman-Bylander form 

At this point we choose an advantageous deviation 
from the original implementation The k-dependent 
non-local potential is arbitrarily complex with respect to 
k. Previously, Shirley expanded the entire operator in the 
optimal basis on a coarse grid in k-space and then inter- 
polated between these values using a polynomial expan- 
sion. This was probably the most complex component 
of the original approach with respect to implementation. 
Details were provided for a quartic interpolation based 
on a 5 X 5 X 5 grid, however, the generalization to differ- 
ent grids and the relative importance of such effort was 
left to the judgement of the reader for specific examples. 
Storage of the ultimate parametrization of the non-local 
potential is proportional to the grid size and the square 
of the number of basis functions. In this work, we take 
a slightly different approach, which we will show to be 
more compact in many cases and more powerful in terms 
of deriving spectroscopic information. 

We assume a generalized, separable, Kleinmann- 
Bylander form for the non-local potential 

V""^ = ^|/3a)^aa'(/3a'| (1) 

AA' 

This is most reminiscent of Vanderbilt's ultrasoft pseu- 
dopotentials [l^. Note that for norm-conserving pseu- 
dopotentials -Daa' is diagonal. The composite index 
A = (J, n, Z,to) refers to the site / of a particular ion 
and its associated atomic quantum numbers. Expanding 
the k-dependent version of this operator in the optimal 
basis, we find that the k-dependence is limited to the 
projectors 

- ^/3A.(k)*i?AA'/9A',(k), (2) 
AA' 

where 

/3A.(k) = (/3A|e'''■'■|B^). 

So, we may consider interpolating only the projec- 
tor matrix elements on a grid in k-space. The k- 
dependent projectors are quite efficiently evaluated by 
one-dimensional Fourier transform of their radial com- 
ponent and should be obtainable from the original elec- 
tronic structure code. In order to make our implemen- 
tation general, we employ three-dimensional B-spHne in- 
terpolation [ia, [m on a uniform ni x 71-2 x 77,3 grid of 



k-points in crystal coordinates, that is, chosen from the 
unit cube, [0,1]'^. Requests for evaluations of the non- 
local projectors at k-points outside the unit cube assume 
periodicity in k-space. In general, we use a larger k-point 
grid to interpolate the projectors than we use for gener- 
ating the optimal basis. A good rule of thumb is to use 
a grid at least twice as dense. We note that for systems 
with c?-electrons we have used more dense grids. We also 
avoid spurious interpolation by limiting the order of the 
B-splines to be equal to the number of grid points in each 
dimension. Note that for an ni x 712 x 713 grid, each di- 
mension is actually expanded by one to include the edges 
of the unit cell. 

Note that by interpolating the projectors, one must 
construct the full non-local potential for each k by matrix 
multipHcation. This apparent additional cost is not that 
great, considering that for norm-conserving pseudopoten- 
tials D\\i is diagonal, and even for ultrasoft pseudopo- 
tentials D\\i is block-diagonal. This specific choice for 
interpolation can reduce the required storage for the non- 
local potential by the ratio of the number of projectors to 
the number of basis functions, which for many applica- 
tions is a reduction on the order of hundreds, given that 
there may be typically 100 basis functions per atom, but 
likely less than 10 projectors. 



B. Extension to ultrasoft pseudopotentials 

The extension of the original approach to ultrasoft 
pseudopotentials is now trivial. Ultrasoft pseudopoten- 
tials relax the norm-conservation condition, by introduc- 
ing a correction derived from atomic all-electron and 
pseudo waves: 

Qa.a'W = VA(r)*VA'(r)-0A(r)*</'A'(r) 



gA,A' = (V^aIV^A') - (0a|0A') 

where -0, 4> refer to all-electron and pseudo waves re- 
spectively. This correction appears in the coefficients 
that define the non-local potential: 

A, A' = i^A.A' + ^A,A' 7 

where the first term is a constant for each atomic 
species, while the second term involves an integral 
of QA,A'(r) over the density dependent Hartree plus 
exchange-correlation potential • 

Furthermore, we must remember that the use of ultra- 
soft pseudopotentials introduces a generaHzed orthonor- 
mality condition on the eigensolutions 

(nk|S'|TOk) = Snm, 



which impHes that their periodic components are solu- 
tions to the following generaHzed eigen-problem: 

[F(k)-e„k5(k)]|u„k) = 0. 

The k-dependent overlap matrix in the optimal basis 
is defined to be 

S,,{k) = ^/3A.(k)*QA,A'/3(k), 

A, A' 

and the projectors matrix elements /3Ai(k) are identical 
with those used in the non-local potential. Note that this 
introduces a large storage and computational saving: we 
need only interpolate the projector matrix elements once, 
and then we can construct both the non-local potential 
and overlap matrix by multiplication. 

In summary, the generalization to ultrasoft pseudopo- 
tentials has the following additional requirements for con- 
struction of the k-dependent Hamiltonian: (1) access to 
the self-consistent coefficients D^ff at the end of the SCF 
calculation and (2) access to the Qx.y coefficients avail- 
able in the pseudopotential definitions. Subsequently, at 
any k-point, one must find the solution of a generalized 
eigen-problem with both i/(k) and S'(k) constructed by 
interpolation. 

C. Advantages for spectroscopy 

One of the most common matrix elements used for op- 
tical spectroscopy is that of the velocity operator. How- 
ever, for non-local Hamiltonians, its evaluation is non- 
trivial, involving a commutator of t he p osition operator 
and the non-local potential [13, [13, H HI, [H ■ Several 
approaches have been introduced to include or overcome 
this complication. For instance, rather than evaluating 
matrix elements of velocity in the transverse gauge, one 
can equivalently evaluate plane-wave matrix elements in 
the longitudinal gauge, exploiting the Heisenberg equa- 
tion of motion: 

(nk|q-v|mk) = lim i^^^^^^±3_i!^(„k|e-*'i--|mk + q>, 

q^Q q 

where q is the polarization of the incident electric field. 
Substitution of the commutator [e~** '', H] in this expres- 
sion leads to the following identity: 



(nk|q-v|TOk) = lim(M„k| l^mk+q) 

q-*0 q 
= (Wrik q Umk)- 

dk 

One advantage of our current implementation is the 
ability to evaluate derivatives of ^^(k) with respect to 
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Figure 3: The atomic structure of 7-brass (CusZng) 



k. The kinetic energy operator is readily differentiable, 
while the B-spline routines we adopted also include effi- 
cient evaluation of derivatives [2l|. Therefore, very little 
extra work is required to access optical transition ampli- 
tudes. Beyond first derivatives, one can see that accurate 
effective masses are easily attainable through the second 
derivatives and would be extremely efficient for large pe- 
riodic nanostructures, obviating the need for numerical 
differentiation based on multiple costly k-point calcula- 
tions. 



V. APPLICATIONS 

Previous work has shown the efficacy of Shirley inter- 
polation in reproducing the band structure of crystalline 
soHds. Instead, we choose to focus on supercells, where 
the efficiency of Shirley's approach is clearly competitive 
with the usual plane- wave calculations, while retaining 
the accuracy of self-consistent results. 

A. Large systems using only the F-point 

We choose 7-brass (CusZng) as an example of a crys- 
talline solid with a large unit cell containing 26 atoms 
(Fig. [3]). We generate the self-consistent field using a 
shifted 4x4x4 k-point grid within DFT using ultrasoft 
pseudopotentials for Cu and Zn, a plane-wave cut-off of 
25 Ry and a charge density cut-off of 200 Ry. The basis 
is built using 200 input states calculated at the F-point, 
which are then expanded to include the seven images of 
the F-point at the corners of the reciprocal space unit 
cube. The optimal basis is obtained by diagonalizing the 
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c 


Npw 


Tpw 


Ms 


Ts 






s 




s 


10 


21993 


191.56 


244 


0.063 


15 


32971 


220.48 


255 


0.056 


20 


43975 


234.88 


279 


0.056 



Table I: Timing information per k-point calculation for 
graphene supercells of varying planar separation c, when us- 
ing a plane-wave DFT code Tpw with a basis of Npw plane- 
waves, and when using Shirley interpolation Ts, with Ms 
basis functions. 



overlap matrix and truncating to 1095 functions from 
a possible 1600, corresponding to ~ 42 basis functions 
per atom. The non-local potential is interpolated on a 
3x3x3 grid. The non-dispersive d-bands of Cu and Zn 
require this level sampling in order to accurately repro- 
duce the non-local potential. The resulting band struc- 
ture is shown in FigureSl Comparison with DFT calcula- 
tions throughout the Brillouin zone indicates remarkable 
accuracy (root-mean-square deviation of 2 meV) for the 
interpolated band structure, with all bands of all charac- 
ter {s,p,d) reproduced to the same degree. We can effi- 
ciently refine a non-self-consistent estimate of the Fermi- 
level using Shirley interpolation, and we find it to be 
shifted from that of our initial self-consistent-field plane- 
wave calculation by 0.3 eV. This is a good illustration 
of the importance of detailed k-point sampling in met- 
als. Shirley interpolation provides an efficient route to 
obtaining a more accurate estimate of the Fermi level 
for a given self-consistent field and indicates the possibil- 
ity of using this interpolation scheme to efficiently refine 
self-consistent-field calculations for metallic systems. 



B. Efficiency in vacuo 

When using plane- waves in supercell simulations of re- 
duced dimensional systems, the inclusion of large vac- 
uum regions comes at a significant computational cost. 
The use of Shirley interpolation can reduce this cost dra- 
matically. We use graphene as an example, where we 
simulate this two-dimensional sheet of carbon atoms in 
a three-dimensional supercell with a large separation be- 
tween periodic images defined by the c-axis. We notice 
that increasing c results in a large increase in the number 
of plane-waves in this dimension, but has only a small 
impact on the number of optimal basis functions used 
to construct the Shirley Hamiltonian. Table U clearly 
illustrates the efficiency of the Shirley approach for k- 
point sampHng. In this small example we see speed-ups 
of greater than 3000. Furthermore, the increase in com- 
putational effort that we expect when adding to the vac- 
uum spacing is practically absent from the interpolated 
case where the number of basis functions increase only 
slightly, rather than the linear increase for plane-waves. 
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Figure 4: The electron band structure of 7-brass (CusZng) 
generated using plane-wave DFT calculations (red dots) and 
using Shirley interpolation based solely on DFT wave func- 
tions generated at the F-point (black lines). The root-mean- 
square deviation of the Shirley interpolated values is 2 meV. 
The black (red) horizontal line indicates the Fermi level re- 
sulting from a 4 X 4 X 4 k-point plane-wave DFT calculation 
(10 X 10 X 10 k-point Shirley interpolation). 



VI. COMPARISON WITH WANNIER 
INTERPOLATION 

In the last decade, the use of maximally localized Wan- 
nier functions (MLWFs) has emerged as an extremely 
efficient and physically appealing route to interpolat- 
ing electronic band structure and deriving useful tight- 
bindin g p arameters from first-principles Hamiltonians 
0, [T3. l27l|. The approach generates Wannier functions 
within a gauge which minimizes their spatial extent. In 
this sense, one constructs a set of basis functions locaHzed 
in real-space, with one MLWF per band. The procedure 
is similar to that outlined in Figure [H with the differ- 
ences lying in an additional minimization of the spread 
of the orthogonaHzed basis functions and no requirement 
to construct a k-dependent Hamiltonian explicitly - this 
is obtained rather by Fourier interpolation. For systems 
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Figure 5: Details of the band structure of 7-brass in three 
energy regions: (a) the Zn d-bands; (b) the Cu d-bands; and 
(c) high in the conduction bands. DFT plane-wave calcula- 
tions are indicated by red dots. Shirley interpolated bands are 
shown as black lines. Together with Fig. |4]this indicates the 
ability of Shirley interpolation to accurately reproduce the 
bands of metallic states of varying character with minimal 
effort (F-point calculations only). 



with bands which do not possess an intrinsically local 
character (sp-bands in metals, for example) the Wannier- 
ization procedure has some intrinsic difficulties related to 
(i) providing a relatively large number of k-points to en- 
able significant localization of the functions and (ii) dis- 
entangling such dispersive bands from manifolds of dif- 
ferent character which they may easily cross due to their 
large dispersion. The latter problem was solved by Souza, 
Marzari, and Vanderbilt pjl, while the former remains 
an intrinsic limitation imposed by the physical properties 
of the system under study. It is particularly problematic 
in spectroscopic studies where large numbers of unoccu- 
pied bands (which are in general dispersive) are needed. 

In contrast with Wannier functions, the optimal ba- 
sis functions used in Shirley interpolation have no con- 
straint on their locaHzation. They are simply the result 
of a diagonalization of the overlap matrix for the entire 
set of input periodic functions. In this sense, generat- 



ing the optimal basis functions is quite automatic and 
does not suffer from issues common to most multidimen- 
sional minimization methods, such as trapping in local 
minima or sensitivity to initial conditions. The resulting 
functions can in fact be quite delocalized and, in gen- 
eral, we have not paid much attention to their spatial 
dependence, given that we do not try to exploit it in any- 
way. For instance, one could not hope to extract tight- 
binding parameters from a set of basis functions which 
are infinitely extended. Figure [6] shows a small number of 
the optimal basis functions derived for fee Cu. They are 
clearly delocalized, and what look like simple functions 
for the larger eigenvalues of the overlap matrix become 
increasingly complex for smaller eigenvalues due to the 
requirements of orthonormality. 

Shirley interpolation is particularly suited to explo- 
ration of metallic band structure, due to the robust au- 
tomatic nature of generating the basis and the obviation 
of disentanghng procedures. Furthermore, one can gen- 
erate the band structure with very few initial k-points. 
In fact, we have already seen that for large supercells the 
F-point is sufficient to generate bands which accurately 
reproduce DFT calculations. Wannier interpolation re- 
quires more k-points to generate accurate band structure 
for bands which do not have an intrinsically localized 
character. For small (monatomic) primitive cells, this 
may not be problematic, but for larger supercells, where 
k-point sampling may still be necessary, then there are 
clear advantages to using Shirley interpolation. 

Finally, it is worth mentioning that a combination of 
Shirley interpolation with the Wannierization procedure 
may be particularly effective for systems with intrinsic 
electron delocalization. Provided that one can generate 
a converged self-consistent charge density, one might use 
Shirley interpolation to efficiently generate solutions to 
the Kohn-Sham equations at as many k-points as desired 
and from these construct the necessary overlap matrix 
elements to begin the Wannierization procedure. This 
may prove particularly advantageous for the interpola- 
tion of metallic or high-energy unoccupied bands in large 
systems, such as metallic alloys, conducting polymers, 
etc. 



VII. FUTURE APPLICATIONS 

The particular advantages of reducing the dimensions 
of a k-dependent Hamiltonian via an optimal basis are 
clear for explorations of band structure and spectroscopy. 
In this section, we outline further possibilities for im- 
proved algorithms or improved scaHng in both DFT and 
beyond-DFT approaches. 

Some self-consistent field calculations rely quite heav- 
ily on numerically converged Brillouin zone integrations. 
For example, an accurate determination of the Fermi- 
level is vital to an accurate estimation of the charge den- 
sity in metaUic systems, and charge transfer at metaUic 
surfaces. For large system sizes, these calculations can 
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Figure 6: A subset of the optimal basis functions of fee Cu, determined using 8 input k-points (F plus seven images) and 
displayed in real-space in a 2 x 2 x 2 supercell. Copper atom positions are indicated by copper-colored spheres. Ouside (inside) 
of basis function isosurfaces indicated in purple (green). Numbers in parentheses indicate the ordering in terms of overlap 
matrix eigenvalue magnitude corresponding to the following coverage of the entire space: (1) 7.6%; (2) 7.3%; (3) 6.4%; (4) 
6.0%; (5) 5.7%; (6) 5.4%; (18) 1.8%; (36) 0.78%. 



prove prohibitively expensive, since the overall cost of 
the calculation scales like NkN^, where Nk is the num- 
ber of k-points and N is the number of basis functions. 
Even though one can assume that larger simulation cells 
reduce the number of required k-points for numerical con- 
vergence, this may not be a sufficient to reduce the overall 
cost significantly. For example, doubling the system size 
would lead to scaling of {Nk/2){2N)^ = 2'^{NkN^), which 
is disheartening if Nk/2 > 1. Since we have seen now that 
for large systems one can quite easily generate accurate 
band energies and states throughout the Brillouin zone, 
one could in principle iteratively calculate just the zone- 
center electronic structure, while using interpolation to 
converge the Fermi-level and self-consistent charge den- 
sity upon which the Kohn-Sham Hamiltonian is based. 
This would reduce the overall scaling, removing the lin- 
ear dependence on Nk at the expense of an increase in the 
overall prefactor associated with generating the optimal 
basis and k-dependent Hamiltonian. 

For calculating excited state properties from fi rst p rin- 
ciples, the combination of the GW approximation [2a and 
Bethe-Salpeter Equation (BSE)f2^, 's^l has emerged as 
an accurate and efficient approach when applied to crys- 
talline solids, molecules, nanostructures, and surfaces. 
This approach is computationally demanding (scaHng at 
least as N"^) and relies heavily on access to detailed Bril- 
louin zone sampling of the calculated electronic structure. 
For periodic systems, a very fine sampling of k-space is 
required to obtain converged BSE solutions, and inter- 
polation procedures have already been applied to enable 
more efficient calculations. In fact, Shirley has already 



used this approach to deal with optical and x-ray exci- 
tations in solids fsl, The main bottle-neck in such 
calculations comes from the need to access the dielectric 
matrix at many k-points. We hope to apply the Shirley 
interpolation scheme to improve the scaling of such cal- 
culations by representing the dielectric matrix within the 
Shirley basis. This will be the subject of future work. 



VIII. CONCLUSIONS 



We have presented a new implementation of the Shirley 
interpolation method. The advances in our approach in- 
clude: (1) a reduction in the number of input electronic 
structure calculations required to construct the optimal 
basis; (2) the ability to interpolate over the entire Bril- 
louin zone using just the zone-center as input for systems 
with large unit cells; and (3) a generalization of the non- 
local potential which reduces storage requirements, per- 
mits the use of both norm-conserving and ultrasoft pseu- 
dopotentials. We provide appHcations of this method 
to sodium, 7-brass, graphene, and copper which illus- 
trate its generality and robustness, particularly in treat- 
ing metals. In this regard, it is competitive with exist- 
ing interpolation schemes based on maximally-localized 
Wannier functions. 
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